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ABSTRACT 

We propose a bilinear sampling algorithm in Green's function 
Monte Carlo for expectation values of operators that do not commute 
with the Hamiltonian and for differences between eigenvalues of dif- 
ferent Hamiltonians. The integral representations of the Schroedinger 
equations are transformed into two equations whose solution has the 
form i(j a (x)t(x, y)tpb(y), where ip a and ^ are the wavefunctions for the 
two related systems and t(x,y) is a kernel chosen to couple x and y. 
The Monte Carlo process, with random walkers on the enlarged con- 
figuration space x (g> y, solves these equations by generating densities 
whose asymptotic form is the above bilinear distribution. With such 
a distribution, exact Monte Carlo estimators can be obtained for the 
expectation values of quantum operators and for energy differences. We 
present results of these methods applied to several test problems, in- 
cluding a model integral equation, and the hydrogen atom. 

Key words: algorithm; bilinear sampling; energy difference; Green's 
function Monte Carlo; quantum expectations; random walk. 
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I. INTRODUCTION 



Quantum Monte Carlo (QMC) methods have found their way into increas- 
ingly many applications in the study of many-body systems. Among them, the 
Green's function Monte CarloW'^'^'M (GFMC) has proved a very powerful way 
to solve the Schroedinger equation in many dimensions. Based on an iterative 
stochastic process, asymptotically it yields a Monte Carlo (MC) representation 
of the ground state wavefunction of the system. For quantum operators whose 
eigenfunctions are the wavefunctions, such as the Hamiltonian itself, only one 
MC sample is needed together with an analytical trial wavefunction to evaluate 
exactly their expectation values. Therefore, the GFMC method is capable of 
calculating ground state energies and related quantities in a relatively straight- 
forward fashion. Indeed numerous highly accurate and effective calculations have 
been performed along these lines for many very different systems. 

Due to the random and discrete nature of MC, however, difficulties arise in 
calculating expectations which involve two wavefunctions. One obvious exam- 
ple is the calculation of expectation values of quantum operators that do not 
commute with the Hamiltonian. One way of calculating expectations with re- 
spect to if) 2 is to generate enough independent samples from ip so that, with 
substantial probability, two lie within range of a known Green's functional . For 
a many-body system, this requires two very large samples. Another example is 
the computation of small energy differences between two systems, described by 
similar Hamiltonians whose ground state eigenvalues differ slightly. The statis- 
tical variance associated with the mean in energies computed independently can 
be comparable to the energy difference, causing the signal in the calculation of 
the difference to be lost in noise. 

Attempts have been made to overcome these difficulties, leading to success 
in specific cases. But the proposed approaches all present certain limitations, 
making such computations either inexact or much more complicated than the 
energy calculation. Among these, the "extrapolation" methodM'^ is the most 
straightforward and has been widely applied to obtain expectation values for 
quantum operators that do not commute with the Hamiltonian. It gives a biased 
estimate to the expectation value, and the bias is often hard to assess. A trial 
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wavefunction is required which accurately describes the desired property of the 
system. Thus for systems that are not well understood, this method will always 
be uncertain. Other approaches^ ' ^ ' ^ ' ^ for expectation values keep track 
of decedents of walkers or employ "side walks" . These methods are in principle 
asymptotically exact and they have been successfully applied to certain problems 
to obtain very accurate results for expectation values. But they are often techni- 
cally quite delicate and asymptotically unstable in the sense that increasing the 
length of the side walks so as to decrease the bias leads to a reduction of signal 
to noise ratio. In the limit of infinite side walks, the ratio is zero. Furthermore, 
their efficiency or even success largely depends also on the quality of the guiding 
wavefunction, i.e., a priori knowledge of the true wavefunction. For the problem 
of energy differences, it is sometimes possible to carefully arrange to correlate 
the two random walks such that the errors in energies would largely cancel^ 12 !. 
But unfortunately this is also limited in its applicability. 

Underlying the difficulties of these calculations is the fact that they require 
highly correlated configurations in the MC process representing two related func- 
tions. We propose in this paper a bilinear sampling algorithm in the GFMC 
framework that achieves such a correlation naturally within a single random 
walk. We show that it is possible to apply this to calculate expectation values 
of quantum operators that do not commute with the Hamiltonian as well as to 
compute energy differences. As a new and very different approach, it results 
from the integral representation of the Schroedinger equations used by standard 
GFMC calculations for energies. Instead of sampling a function linear in the 
unknown wavefunctions, we transform the Schroedinger equations into a pair of 
integral equations whose solutions are bilinear in the wavefunctions. The ran- 
dom walk based on the new pair of equations directly samples the ground state of 
such solutions, namely, ip aQ (x)t(x,y)i^bo(y), where a and b label the two related 
systems and ifj aQ and t/^q are correspondingly their ground state wavefunctions. 
(The labels a and b can be the same, in which case the wavefunctions describe 
the same system.) The function t(x, y) is chosen to couple appropriately the 
configurations x and y. The object of the random walk is an ensemble of pairs 
of configurations (x, y) rather than individual configurations x. 
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As we shall show in Section II, the bilinear sampling method yields asymp- 
totic distributions of configurations (x,y) which can provide exact MC estima- 
tors for the function tjj a0 (x)ij)bo( x )- When a = 6, this is simply the square of 
the ground state wavefunction for a system, and we thus have an exact way 
to compute ground state expectation values of quantum operators. If a and b 
are different, it permits a direct calculation of the energy difference as AE = 
< *p aQ \ AHfybQ > I < ipaol^bo >> where AH = H a — is the difference in the 
Hamiltonians for the two systems. 

As tests, we have applied this new algorithm to a model problem to com- 
pute various moments and also to the ground states of the hydrogen atom and 
related systems. The former is based on an integral equation composed entirely 
of Gaussians. Because of its simplicity and flexibility, it provides a very trans- 
parent picture of various aspects of the problem and enables us to study the 
algorithm from many different angles. In the latter, we calculate expectation 
values of various operators for the ground state of H and also energy differences 
between the ground states of H and similar systems with different potentials. 
This provides tests of the necessary ingredients for the application of the new 
algorithm to larger systems of the same class. 

In Section II, we outline the formalism of the bilinear sampling method 
and its numerical implementation. Then the application to two model systems is 
described in Section III. In Section IV, we give a discussion of the method and its 
further possibilities. Finally in the Appendix the sampling techniques involved 
in the hydrogen calculation are developed. These are in fact general techniques 
necessary for a class of such calculations for atoms and molecules. 

II. FORMALISM AND THE RANDOM WALK 

As we shall discuss in Section IV, the bilinear sampling approach is by no 
means limited to a certain type of Green's function. But for simplicity in presen- 
tation, we use the original form of GFMC^. In atomic units, the Schroedinger 
equation for a many-electron system can be written as 

[-±V 2 + V(xM(x)=E^(x), (1) 
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where a; is a 3M dimensional vector denoting the coordinates of all M electrons 
in 3D space, and the energy is negative for any bound state. We will use the 
Green's function g(x,z) for the operator (— -|V 2 — E) in 3M dimensions. The 
energy E is unknown, but it can be either scaled away for Coulomb systems! 13 ] 
or obtained iteratively, which usually converges very fast (see Section III. 2. 2). 
We can then transform the Schroedinger equation (1) into the following integral 
equation: 

ip(x) = A / g{x, z)w(z)ip(z)dz. (2) 



The Green's function g here has an analytically known form and w = —V. Eq 
(2) has a set of solutions ip with different eigenvalues A, the lowest of which is 1. 
With an arbitrary initial function having non-zero overlap with the ground state 
wavefunction, this equation can be iterated to yield asymptotically the solution 
corresponding to the lowest A, or the ground state ipo. This forms the basis of the 
GFMC method. In practice, the process is carried out with the wavefunction in 
each iteration represented by a generation of individual configurations, or random 
walkers. The walkers live in configuration space and they move from one point 
z in this space to another x according to the probability distribution function 
g(x,z). The function w in Eq (2) is treated as a multiplicative weight or as a 
source of branching of walkers. 

We consider two related systems described by Eq (2): 

ipa{x) = A a J g a (x,u)w a (u) ip a (u)du (3a) 
tp b (y) = \ b J g b (y,v)w b (v) ip b (v)dv, (36) 

where again the subscripts a and b denote the two systems respectively. In order 
to obtain an asymptotic solution bilinear (rather than linear) in the wavefunc- 
tions, we use a coupling function t(x, y) and multiply each equation in Eqs (3) 
by t(x,y) and the wavefunction for the other system to arrive at the following 
pair of equations 



,(x)t(x,y)ip b (y) = X a J t X, ^ a y j ,U w a{u) ip a (u)t(u,y)il; b (y)du (4a) 
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and 

ip a (x)t(x,y)ijj b (y) = \ b I t -^j^^y^- w h {v) ^ a {x)t{x,v)^ b [y)dv. (46) 

J tyx, v) 

These hold for any positive coupling function t in principle and their solution 
now has the form ip a (x)t(x,y)ip b (y). Eqs (4) are completely defined once the 
kernel t is chosen. In this paper, we assume that t is symmetric in x and y solely 
for simplicity. Eqs (4) can be rewritten so that they are more transparent for a 
random walk interpretation: 

®(x,y) = X a J ' T a (x,y\u,v)N a (u,v) $(u,v)dudv (5a) 

and 

$>(x,y) = X b J T b (y, x\v, u)N b (v, u) &(u,v)dudv, (56) 

where $>(x, y) = ip a (x)t(x,y)ipb(y) is the bilinear solution we seek. The kernel 
T s is a normalized probability distribution function of x and y conditional on u 
and v defined as 



F s (x, y\u, v) oc [g s (x, u)t(x, v)]S(y - v), (6) 

where s is either a or b. The multiplicative factor iV s is 

N s (u,v) = w s (u) J g s (x,u)t(x,v)dx/t(u,v). (7) 

We note that neither T s nor iV s is symmetric in its variables. 

To solve Eqs (5) by Monte Carlo for the ground states in the bilinear form 
$o = i ; ao(x)t(x, y)iibo(y), we introduce random walks on the enlarged configura- 
tion space x®y. In other words, each walker now consists of two configurations, 
namely x and y, which sample a joint distribution function 3>(x, y). The itera- 
tion, or the random walk, is carried out by moving walkers according to either 
one of the two equations in each step. The initial distribution of walkers can be 
arbitrary or generated from a Metroplis sampling of $ with the -i/Vs replaced by 
trial wavefunctions for the ground states. The bilinear sampling technique en- 
ables us to treat explicitly functions bilinear in the wavefunctions and maintain 
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highly correlated configurations. But since Eqs (5) are a direct transformation of 
Eqs (3), the convergence property of the original linear integral equations to the 
corresponding ground states simply carries through. In the random walk process, 
both equations need to be applied about equally often in order to ensure conver- 
gence and assure efficiency. In practice, this can be accomplished by using the 
two equations alternately, or randomly with equal probability. Once an equation 
has been selected in a step, the multiplicative factor (either N a (u, v) or Nb(v, u) 
as determined by the equation) is constructed for each walker. New walkers for 
the next generation are then produced by: 1) choosing a walker (u, v) from the 
old population either with probability proportional to their multiplicative factors 
or by branching (depending on whether the size of the population is fixed or not); 
and 2) sampling a new walker from the parent walker according to the kernel T s 
for that equation. We note from the form of F s that in fact in 2) only one new 
configuration (x or y) is selected and the complementary configuration in the old 
walker is simply carried along. 

Much freedom still remains about the choice of the coupling function t. If we 
set t to a constant, bilinear sampling would reduce to generating two independent 
sets of configurations from two standard GFMC runs, from which the overlap 
would be difficult to extract in high dimensional systems. The kernel t must 
generate configurations that are close together. The requirement of this is very 
clear in the current problem. When t is the same as one of the Green's functions, 
say, g a , the bilinear sampling method yields a density ^ aQ (x)g a (x, y)i^bo(y)- From 
the original integral equation, Eq (3a), we have 



As mentioned in Section I, the product on the left-hand side of Eq (8) is essential 
to the expressions for the expectation values we seek. In the asymptotic regime 
of the MC process, the distribution of pairs of configurations represents 3>o(x, 2/), 
as a sum of delta functions in x and y: 







k 
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where {x k ,yk} is the collection of these pairs (random walkers) labeled by k. 
This combined with Eq (8) gives, 

^ao(y)^b (y) oc ^2w(x k )5(y -y k ). (10) 

k 

In other words, w a (xk) is an MC estimator for iftaoiyktyboiUk) ■ 

With the direct sampling of the product of two wavefunctions given by Eq 

(10) , we can easily obtain the desired expectation values of quantum operators 
exactly. If a = 6, the ground state expectation value of any multiplicative quan- 
tum operator O is 

<Cy - ! ^*o(y)°(y)^ao(y) d y _ E fc 0(^)^0^) 

I ^a {y)ipa {y)dy Y,k w ^k) 

When a ^ 6, we wish to calculate the ground state energy difference between 
two systems a and b described by different but related Hamiltonians H a and H b . 
Their ground state wavefunctions are given by ifj a0 and ip bQ , which are assumed 
to be non-orthogonal to each other. If H a — H b = V a — V b , since 

AE = J ^ao(y)H a ^ bo (y)dy _ J ip aQ {y)H h ip hQ (y)dy 

J ^a (y)^b (y)dy J ipao(y)^b (y)dy 

_ I ^ao(y)[Hg - H b ]if; b0 (y) ^ 
I tpa (y)ipb (y)dy 

we have 

AE = T /k w a( x k)[Vg(yk) - Vbjyk)] ^ 

Of course, in both cases we can use Eq (3b) in Eq (8) to obtain w b (yk) 

as an estimator for ipao(xk)^bo( x k), which can provide expressions similar to Eqs 

(11) and (13). They can be combined with the equations shown here for better 
statistics. 

Implicit above is the assumption that it is possible to evaluate N as well 
as to sample V. This is, of course, not universally valid. Fortunately, in the 
current GFMC approach, the Green's function for (— -|V 2 — E) for any system 
can be written as a superposition of Gaussians and it is straightforward to sample 
product of two such functions. The integrals in N can also be easily obtained. 
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The assumption is also true for other classes of Green's functions for various 
systems of interest. For discrete systems such as certain quantum spin systems, 
it will be even less challenging in principle. 

III. APPLICATIONS 

III.l. A Model Problem 

In this part, we apply the bilinear sampling method to a model problem^ 14 ! 
described by an integral equation of exactly the same form as the general many- 
dimensional equation given by Eq (2). We shall test the capability of the al- 
gorithm in calculating expectation values by evaluating various moments of the 
"ground state" distribution. The labels a and b (thus s) may be omitted in this 
section. Without altering the notation, we simply redefine as follows 

g(x, z) = exp[— a(x — z) 2 ] 

and 

W{Z) = V^T eXP( "4o^2 22) ' 
where a can be any real number greater than 0.5. It can be easily verified that 
under this new definition, Eq (2) has a set of solutions which are product of 
Gaussians and polynomials. Among them the ground state, the solution corre- 
sponding to the lowest eigenvalue A = 1, is 

if> (x) = exp(-^a; 2 ). 

We choose the coupling function t to be also a Gaussian, 

t(x,y) = exp[-P(x - y) 2 ], 

and we shall solve Eqs (5) by MC for the joint distribution $ of the ground 
state. 

With every term in Gaussian form, this problem is easy to study analytically 
from every aspect. The parameters a and (3 are completely at our disposal and 
they can be varied to provide insight into the behavior of the algorithm under 
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very different conditions of the Green's function and coupling. Furthermore, if 
necessary, the iteration can also be carried out directly without doing MC. That 
is, assuming a general solution exp(— a n x 2 + b n xy — c n y 2 ) for the n th iterate 
of Eqs (5), we can determine a n +i, and c n+ i at every stage. This can 

be used to generate numerical trajectories so as to observe the convergence. 
The product g(x, u)t(x, v) in kernel F(x,y\u,v) can be easily transformed into 
a single Gaussian in the unknown x by completing the square in the exponents. 
Not surprisingly, then, sampling V amounts to sampling a Gaussian. 

The MC process generates random walkers {xk, yu} representing the distri- 
bution if;o(x)t(x,y)^o(y) for the ground state. Since 

J ij>o(x)t(x, y)ipo(y)exp(-j^—^y 2 )dy oc ip%(x) = exp(-x 2 ), 

the MC result is adequate to determine completely the function exp(— x 2 ). For 
example, various moments can be computed in a similar fashion as in Eq (11) 
and compared with the exact results. Also a histogram can be easily made for 
this function in one or two dimensions. The second moment < x 2 > is always 
computed in our tests. We have carried out calculations in nine dimensions and 
the result is satisfactory. For the present study, however, one dimension suffices 
in revealing the characteristics of the algorithm. A wide range of values have 
been used for both a and (3. Table I presents our results for the second moment 
in one dimension. We see that all the computed answers agree well with the 
exact result, namely 0.5. 

One important issue is how the algorithm performs when a and (3 are large, 
since that is the case when both kernels are very sharply peaked. That the 
Green's function is sharp implies the step size of the random walk is small in 
configuration space, which more closely resembles the situation in higher dimen- 
sions resulting from actual many-dimensional problems and also the situation in 
which a QMC is generated by a diffusion process. The coupling is consequently 
also very sharp, as indicated in Section II. In fact in many cases such as the 
hydrogen atom below, the coupling is simply the Green's function. Under this 
circumstance, the probability becomes extremely small to have two independent 
and random configurations appear close in configuration space. Therefore the bi- 
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linear sampling method must effectively couple two sets of configurations without 
distorting the distribution. It is reasonable to expect less efficiency as a and (3 
increase, but the sampled distribution must nevertheless be correct. From Table 
I, we see this is indeed the case. 

We keep a constant number of walkers in our simulations. This results in a 
bias, since generations with high multiplicity contribute less than they should and 
vice versa. Branching can be introduced instead of strict population control to 
avoid this. It is also possible to correct for such a bias by carrying weights^ 15 !. The 
sum of multiplicative factors over all walkers in a generation labeled by n, J\f n , 
indicates the total number of walkers the next generation should include. Thus 
with the number of walkers fixed at L for every generation, each generation can be 
assigned a weight formed by a product of Af/L from a certain number of previous 
generations, i.e., W n = Y\aL\ N n -i/L. The number of previous generations to be 
included, m, can be tuned so that it is large enough to remove the bias and yet no 
excessive fluctuation is introduced. For this model problem, with a fairly small 
number of walkers (usually 1000), it requires less than ten generations to correct 
for the bias. It is observed that, without the correction, the bias effect is often 
quite significant. As the kernels become sharper, the bias becomes more and 
more serious. Moreover, it does not seem to always exhibit a clear 1/L behavior 
as in many GFMC calculations. The population control bias can be attributed 
to fluctuation of the multiplicative function N. In bilinear sampling, due to 
the extra function t inserted to couple two points, it is not implausible to have 
relative large variation in N. For instance, in this model problem, N(u, v) oc 
exp[/3 2 (-u — v ) 2 / (a + (3)]w(u). In these calculations caution must be exercised to 
ensure that the result is unbiased, in other words, robust against population size. 

We also mention in passing that even though Eqs (4) and therefore Eqs (5) 
are true for any non-zero t, they are not necessarily always well-behaved in a MC 
calculation. To illustrate this, recall the variance^ 16 ' of the total weights for a 
generation in the random walk is given by f N(u, v) 2 Q(u,v)dudv. But when the 
coupling function t is much sharper than the Green's function g, this expression 
diverges, which implies that the MC sampling would not actually converge to a 
definite answer. For the model problem, it is possible to determine the range 
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of (3 for each a where infinite variance can be expected. We have verified that 
when the parameters are given in this range, the MC answer can disagree with 
the exact one significantly. The above analysis has employed no knowledge of 
the specific form of the kernels and thus is general to bilinear sampling. For the 
purpose of computing the expectation values of quantum operators, however, it 
is always possible to avoid a kernel t in that regime. So this should not pose any 
problem. As an additional probe of the diverging weights, we can monitor the 
fluctuation of population sizes in a calculation. 

III. 2. The Hydrogen Atom 

As another test case, we use the algorithm to study the ground state of the 
hydrogen atom. In the first part, expectation values of operators with respect 
to the ground state are computed. In the second, we study similar systems with 
different potentials and evaluate the energy difference between the ground states 
of the new and original systems. It is assumed the systems are non-relativistic 
with the nucleus fixed. Except for technical details, implementation of the new 
algorithm to address these problems directly follows the formalism developed in 
Section II. Based on Eq (8), the natural choice of t is t = g a , where g a is the 
Green's function for the hydrogen atom. 

In general, the Green's function g for an M-electron system as defined in 
Section II can be written in an integral representation, which is a superposition 
of Gaussians with different widths. There also exists an analytical expression for 
g in terms of polynomials and modified Bessel functions of the second kind so 
that it can be conveniently evaluated. Using these expressions, the integral in N 
can be easily computed for any M and there exists an efficient way to sample x 
from the product of functions in the same class, g a (x,u)gb(x, v) (cf. Appendix). 

III. 2.1 Ground State Expectation Values 

Similarly to III. 1 , we can sample the square of the ground state wavefunction 
and compute directly within a single run expectation values of multiplicative 
operators. The sampling techniques are only a special case of that discussed in 
the Appendix. In Table II, we show computed expectation values of the potential 
energy V, the radial distance x 2 , and the square of the third component of 
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the electron coordinate together with the exact answers. The number of 
walkers is typically 3000 and the bias is not noticeable. Because of the simplicity 
of the sampling process and the direct sampling of the product of wavefunctions, 
the code is fast and the algorithm quite efficient. We see that the agreement 
between the bilinear sampling and the exact results is again excellent. 

III. 2. 2 Energy Difference Calculations 

In this part, we describe tests of the capability of the bilinear sampling 
method to calculate energy differences by considering systems similar to the H 
atom but with different potentials. We will employ bilinear sampling to compute 
the energy differences between the ground states of such systems and the hydro- 
gen atom. The Green's functions corresponding to the two Hamiltonians belong 
to the same class discussed in the Appendix, only with a possible difference be- 
tween multiplicative factors in their arguments that is a function of the energy. 
The difference in the potential terms causes w a to differ from Wf, in Eqs (4). 

The simplest way to obtain such a system is to add a perturbation term to 
the original hydrogen Hamiltonian, i.e., Hb = H a — ^H', where H' is a multiplica- 
tive operator and 7 is a small coefficient. From Eq (13), we can compute AE/j 
as a function of 7, where AE is the energy difference between the ground states 
of the two systems. When 7 = 0, this in fact is exactly the bilinear approach of 
calculating the expectation value of H' in the ground state of H a . Therefore, in 
a sense, all results in III. 2.1 can be viewed as special cases of these calculations. 
As 7 is increased, the result should deviate from the ground state expectation of 
H' and should always give the exact energy difference divided by 7. From the 
MC point of view, it implies the ability to compute the exact energy difference 
for all ranges of the parameter 7. Moreover, the effect of the small perturbation 
is generated with small fluctuation. In fact, the statistical error may decrease 
with 7. 

We consider a perturbation H' = l/\x\. This is just the original Coulomb 
potential and the energy difference is trivially obtained analytically. In Fig 1, 
we plot the computed AE/'j for some values of 7 from bilinear sampling and 
compare them with the exact result. The ground state energy of H is —0.5 in 
atomic units. We see that the agreement is excellent. For instance, in the case 
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of 7 = 0.003 the bilinear sampling method easily yields an energy difference of 
0.0030036(12), which would be extremely challenging, if possible at all, for an 
approach by independent MC calculations. 

We next study a system described by the Hamiltonian Hf, for the so-called 
Hulthen potential 

t r h T/ exp(-p|a;|) 
1 — exp(— p\x\) 

where Vq and p are parameters. This system can be solved exactly! 17 ] and its 
ground state energy is given by Et, = — (2Vo/p — p) 2 /8 (p 2 < 2Vq). Let 
Vq = p. Then this potential behaves like the Coulomb potential at small values 
of \x\ and approaches zero exponentially at large distances. By varying the 
parameter p, we can control how similar this system is to the hydrogen atom 
and the energy difference between their ground states can be calculated from the 
bilinear sampling method and compared with the exact results. In Table III, we 
show the computed and analytical results for AE for some values of p. Again 
extremely accurate values are easily obtained with bilinear sampling for small as 
well as large energy differences. (Each number corresponds to roughly one hour 
on an IBM RS6000 workstation.) 

As mentioned above, the argument of the Green's functions scales with the 
ground state energies and have the form g(k s \x — y\), where g is the Green's 
function for the operator (—V 2 + 1) and k s = y/2\E s0 \ (s = a,b). This does 
not pose any difficulty because we can obtain iteratively and very quickly the 
correct values for the two energies. In fact, it is observed that, quite generally, the 
final result of a standard GFMC calculation in this approach is rather insensitive 
to the initial input of the energy. We have tested the effect of iterations from 
a very poor starting value of the energy Ebo in our calculations with the l/\x\ 
perturbation and indeed, the convergence is very fast. For example, if we use E a0 
as initial value for both Green's functions, then within one run, we can obtain 
the first order correction to E(, Q . Even for large 7 the effect of the uncertainty of 
an initial energy guess becomes unnoticeable in the final answer in one to a few 
more iterations. In Fig 2, we illustrate this by plotting AE/^y as a function of 
iterations for two large values of the coefficient 7 (0.1 and 0.5). 
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The number of walkers in these calculations is typically 2000 and we correct 
for the bias with several previous generations as described in Section III.l. We 
alternate the two equations in Eqs (5). In the bias correction, we calculate the 
averages of Mi separately for even and odd / and normalize each accordingly such 
that these multiplicity factors are kept around unity. 

IV. DISCUSSION 

It is interesting to note that the manipulation of the original integral equa- 
tion (3a) to arrive at Eq (4a) is very similar to an importance sampling 
transformation' 1 ] 'M. In fact, if the coupling kernel t(x, z) is the same as g a (x, z), 
we can use as importance function the correct form, namely, the unknown wave- 
function ip a0 in its integral representation as given by Eq (3a). Since the MC 
process in effect does the integral in this importance function, we only use the 
integrand and also drop the potential term w a and Eq (4a) ensues. This also sug- 
gests that the algorithm should be rather efficient and, explains to some degree 
why no trial wavefunction is needed in the method. 

As a straightforward method that depends little on a priori knowledge, bi- 
linear sampling should find itself useful to different quantum problems as the 
GFMC approach is employed more and more to understand various many-body 
systems. For example, it seems possible to study with this method ground state 
properties of certain quantum spin systems' 11 ]. A generalized version of the ideas 
developed here with quadrilinear sampling provides the possibility of calculating 
transition moments between two quantum states, as opposed to a method using 
side walkst 10 !. 

By sampling the product of wavefunctions, the bilinear approach is more 
promising than dealing with two independent sets of configurations' 5 ]'' 10 ] from 
two GFMC calculations. We need to study more the behavior of the method as g 
and t become very sharp so as to give a general prescription for avoiding the large 
biases or fluctuations we have seen in that limit. In many-electron systems, we 
inevitably will encounter the "sign" problem in quantum MC^'t 13 L We have not 
yet formulated a bilinear sampling algorithm mechanism that also addresses that 
problem in an exact way. But it is possible to use bilinear sampling within the 
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fixed-node^ 3 !'! 18 ! approximation. With this, the method is clearly generalizable 
to many-electron systems in which the random walk is either generated by a 
diffusion approximation or by domain Green's function methods. In either case, 
there are three possible choices for t(x,y). We can use g{x,y) as in the present 
work, ignoring the fact that a small number of estimates will be negative (because 
w(x) = —V(x) is negative.) We can couple the two configurations x and y using 
the kernel by which new points are generated (itself a Green's function either in 
a short time approximation or over some finite domain.) Finally, it is possible, in 
principle to couple the walkers by the full Green's function (unknown in advance, 
but generated by walks that would go from x to y) and modifying the coupling 
recursively using local Green's functions. Investigation of these alternatives will 
be the subject of future research. 

The bilinear sampling algorithm derives directly from the integral represen- 
tation of the many-body Schroedinger equation of the system. It appears to 
be quite natural for the problems involving functions quadratic in ground state 
wavefunctions, and for the calculation of energy differences. 
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APPENDIX: SAMPLING THE PRODUCT OF TWO GREEN'S FUNCTIONS 



In this appendix, we shall complete the technical part associated with the 
sampling of kernel T s in Eq (6) and the evaluation of the multiplicative factor N 
in Eq (7). These are general to Green's functions for the operator (— |V 2 + \E\) 
with any number of particles. As discussed in Section III. 2. 2, we can assume the 
energy is known. 

Again let M be the number of particles in the system we treat. Then the 
Green's functional'! 19 ] as defined above is given by 



g(x, z) = (^? M/2 [ t-^eM-t ~ k -^ L )dt, (Al) 

where k = y/2\E\. As already mentioned, g also has a form which can be used 
to evaluate the Green's function: 

g(x, z) = (^) 3M / 2 K 3M/2 _ 1 (fc|x - z \)/{k\x - z\f M ' 2 ~\ (A2) 

where K m is the modified Bessel function of the second kind. In the bilinear 
sampling process we need to evaluate N s and sample Y s . With t(x, y) = g a (x, y), 
both N s and T s involve the product of two Green's functions of the same class. 
It is necessary to sample such products as well as evaluate integrals based on 
them. 

We consider the product of two Green's functions with k a and k^. Let £ = 
(^6 — ^ 2 )/2 (> 0) and a = (/c 2 + /c 2 )/2. From Eq (Al), by completing squares 
and changing variables, it is straightforward to obtain the following expression: 

9b(x,y)g a (x,z) 

r°° f p (v — z) 2 

oc^ J T(x\y,z,r u r 2 )exp(-^) P~ 3M/2 exp(-ap - ^ - } ) dqdp. (A3) 
The new pair of variables p and q are 

p = Ti + r 2 
q = ri - r 2 



19 



and T(x\y, z, t±, t 2 ) is a normalized probability distribution function of x condi- 
tional on 7-^2 (or p and g), y, and z: 

T( X \y,z,n, T2) = { ^f^e M -Il^ {x -yil±^n (A4) 

47TT1T2 4TiT 2 T 2 + Ti 

Now we describe the actual sampling and integrating of this product of 
Green's functions gb(x, y)g a (x, z) given the positions of the parent walkers y and 
z. Since integration over x in Eq (A3) simply removes T, we have the following 

J g b (x, y)g a (x, z)dx oc ^[k 3 a M - 2 g a (y, z) - k 3 b M - 2 g b (y, z)]. (A5) 

As a special case, in the limit k a = k bl i.e., when the two Green's functions 
correspond to the same system, the above expression reduces to: 



/ 



g a (x,y)g a (x,z)dx oc ^M" 4 ^ 2 " 1 ^ z), (A6) 



where the superscript of g indicates the Green's function is for 3M/2 — 1 di- 
mensions rather than 3M/2 for the original functions. To sample an x from a 
probability distribution function proportional to g b (x,y)g a (x, z), we note that 
for any known pair of y and z, Eq (A3) can be written as 

9b(x,y)g a (x,z) oc J J T(x\p, q)Q(q\p)P(p)dpdq, (A7) 

where T is the same as in Eq (A4) with the variables changed from 7~i, r 2 to p 
and q. Q is a normalized probability distribution function of q conditional on p 
as follows 

Q(q\p) = { £ ex P(~^)/[ ex P(£P) - exp(-fp)], if \q\ < p ^ 
vwifv I q ? otherwise. 

P is a positive function of p on (0, 00) which therefore can be viewed as a prob- 
ability density function 

l-exp(-2£p) _ 3M/2+1 . 2 (2/-^) 2 n / An ^ 
^(P) oc — p / + exp(-k a p — — ). (A9) 

Thus to obtain x according to the probability density given by Eq (A3) or Eq 
(A7), we need to sample a p from Eq (A9), then agon [— p, p] according to (A8) 
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and finally sample the Gaussian in Eq (A4). Sampling of (A9) is elementary t 16 l 
and can be accomplished by, e.g., sampling the exponential distribution in p 
in the last term and then doing rejections. The case k a = kb is once again 
straightforward, since Q becomes uniform on [— p, p] and the first term in P(p) is 
simply 1. We also mention that analogous (though possibly less elegant) methods 
will apply for any coupling kernel that can be written in the form 

f°° (x-v) 2 
t{x,y) = J^ h{t)eM~ J )dt 

for h(t) > 0. 
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Table I. 



Oi 






■ 


result 






0.6 


0.6 : 


5005(5) 


1.0 : 


5001(4) 


3.0 : 


501(6) 


1.0 


0.6 : 


0.5001(2) 


1.0 : 


0.4998(2) 


3.0 : 


0.503(5) 


3.0 


1.0 : 


0.4999(7) 


3.0 


. 0.501(3) 


5.0 : 


0.501(2) 


5.5 


1.0 : 


0.5003(9) 


3.0 


. 0.500(4) 


5.5 : 


0.502(3) 


10.5 


0.6 


: 0.500(1) 


3.0 


. 0.501(3) 


10.0 


: 0.504(5) 



Table Caption 

Table I. Results of the bilinear sampling method applied to the ID model prob- 
lem. Shown is the second moment < x 2 > from the sampled distribution for 
■i(j 2 (x). The exact answer is 0.5. The parameters a and (3 give the sharpness of 
the Gaussian kernels g (the Green's function) and t (the coupling) respectively. 
The statistical errors in the results are in the last digits and are indicated in 
parentheses. 
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Table II. 



item 


< V > 


< Vx 2 > 


< x 2 > 


< 4 > 


bilinear 


-1.001(2) 


1.500(3) 


3.000(9) 


1.002(3) 


exact 


-1.0 


1.5 


3.0 


1.0 



Table Caption 

Table II. Results of the new algorithm applied to the ground state of the hydrogen 
atom together with the exact answers. The items are the expectation values of 
the potential, the radial distance, the second moment, and the ^-component of 
the second moment. All quantities are in atomic units. The statistical errors in 
the MC results are in the last digits and are indicated in parentheses. 
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Table III. 



p 


exact 


bilinear— exact 


0.001 


0.00049987500 


-.00000000004(8) 


0.0125 


0.00623047 


.00000000(1) 


0.03333 


0.01652778 


.00000007(6) 


0.08333 


0.0407986 


-.0000001(4) 


0.4 


0.180000 


-.000004(4) 



Table Caption 

Table III. Energy differences AE = Eb — E a0 from bilinear sampling between 
the ground states of the hydrogen atom and a similar system described by the 
Hulthen potential, compared with exact results. The Hulthen potential is given 
by Vhui(M) = — pexp(— p|x|)/[l — exp(— whereas the Coulomb potential 
in H is V^dxl) = —l/\x\. Atomic units are used. The ground state energy of H 
is E a0 = —0.5 and Eb is higher. The first column is the exact result for AE, 
while the second column gives the error (bilinear— exact) in the bilinear result. 
The statistical errors from MC are again in the last digits and are shown in 
parentheses in the second column. 
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Figure Captions 

Fig. 1. The "effective" energy difference AE/'j between the original and 
perturbed systems for a perturbation ~/V(x) to the ground state of the hydrogen 
atom. Results of the bilinear sampling method are shown with statistical fluc- 
tuations and are compared with the exact result. V(x) is the potential energy 
operator. The bilinear calculations were done with the arbitrarily chosen values 
7 = 0.003,0.005,0.01,0.05,0.1. The first order perturbation result is 1. Atomic 
units are used. 

Fig. 2. Convergence of iterations with the input of the energy Eb of the 
perturbed system in the energy difference calculations for H. The perturbation 
is again — 'y/jarj and large perturbations (7 = 0.1 and 0.5) are chosen in order to 
show a visible convergence process. Values of AE/'j are plotted together with 
statistical errors as a function of the number of iterations. The exact results 
are given as straight lines. In both calculations an initial value is assumed 
for the energy difference. In one iteration, they both give AE/j = 0.9995(9), 
equivalent to the first order perturbation result. The calculation with 7 = 0.1 
(lower curve) requires only one more iteration to converge to the exact value 1.05 
while that with 7 = 0.5 converges in four more iterations to the correct answer 
AE/7 = 1.25. 
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